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Abstract. The time evolution of £-spin reduced density operators is studied for 
a class of Heisenberg-type quantum spin models with long-range interactions. 
In the framework of the quantum Bogoliubov-Born-Green-Kirkwood-Yvon 
(BBGKY) hierarchy, we introduce an unconventional representation, different 
from the usual cluster expansion, which casts the hierarchy into the form of 
a second-order recursion. This structure suggests a scaling of the expansion 
coefficients and the corresponding time scales in powers of N 1 / 2 with the system 
size N, implying a separation of time scales in the large system limit. For 
special parameter values and initial conditions, we can show analytically that 
closing the BBGKY hierarchy by neglecting £-spin correlations does never lead 
to equilibration, but gives rise to quasi-periodic time evolution with at most 1/2 
independent frequencies. Moreover, for the same special parameter values and 
in the large- N limit, we solve the complete recursion relation (the full BBGKY 
hierarchy), observing a superexponential decay to equilibrium in rescaled time 
r = tN' 1 ' 2 . 



1. Introduction 

Equilibration and thermalization in closed quantum systems have recently seen 
renewed interest, triggered by the impressive progress in performing experiments with 
ultracold atoms and ions [T|. In these experiments the coupling to the environment is 
negligible on the accessible timescales, and they can therefore be regarded to a very 
good approximation as closed systems. Typically, equilibration (in a suitable sense) is 
expected to take place in closed quantum systems [2j El E] , but important exceptions 
are known to exist. Failure of equilibration in systems close to integrability has been 
observed experimentally by Kinoshita et al. [3] , and a substantial body of theoretical 
and numerical work has been devoted to this problem over the last years (see [6] for 
a review). 

Studying the long-time dynamics of many-body quantum systems in general is 
a daunting task. For numerical work, one is usually restricted to fairly small system 
sizes, much smaller than most experimental realizations and not anywhere close to the 
macroscopic regime. The situation is different for integrable systems where analytic 
solutions may exist and larger systems may be studied. However, these systems 
are known to show atypical behavior, possibly very different from the non-integrable 
systems one is often interested in [7]. 
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One solvable model (a quantum spin model with Ising-type interactions and 
subjected to an external magnetic field) for which the long-time evolution can be 
studied analytically, was proposed by Emch [5] in 1966 and later studied by Radin 
[9] in more detail. For a certain class of initial conditions, the time evolution of the 
expectation values of operators A of a certain kind can be calculated analytically for 
arbitrary system sizes N, showing non-Markovian relaxation to the thermal average 
in the thermodynamic limit TV — > oo [8]. This model has been extended by one of us 
to the case of long-range interactions, i.e. spin-spin interactions decaying like r~ a with 
the distance r on a d-dimensional lattice, where the exponent a satisfies ^ a ^ d 
|10U11| . As in the short-range case, the expectation values of the operators A relax to 
thermal equilibrium, but they do so on a time scale that diverges with the system size 
N: The larger the system is, the longer it takes to thermalize. In fact, the relaxation 
dynamics becomes so slow that, within a given experimental resolution and for large 
enough N, no deviation from the initial state will be observed. 

A similar phenomenon, going under the name of quasi-stationary states, has 
received considerable attention in the field of classical long-range interacting systems 
|12[ [13] . and in particular in classical gravity [14]. In the classical context, kinetic 
theory, and in particular the Vlasov (or the collisionless Boltzmann) equation, has 
been successfully applied to describe these quasi-stationary states [151 [T6l [l~2l US] [IT] . 
Such a Vlasov description is known to become exact for long-range systems in the 
thermodynamic limit, and is therefore also expected to be a good starting point for 
a description of large but finite systems |18] . The main goal of the present work is 
to study, similarly to the classical case, the long-time persistence of quasi-stationary 
non-equilibrium states and their final relaxation to equilibrium for long-range quantum 
spin systems within the framework of quantum kinetic theory. We can then test the 
predictions of the quantum kinetic theory against the analytic result for the Emch- 
Radin model reported in |1Q[ [TT] , but the scope of such a kinetic theory goes beyond 
this model: On the basis of a tested and well-working kinetic theory, equilibration 
can then be studied for non-integrable generalizations of the Emch-Radin model, or 
non-integrable long-range variants of the Heisenberg model for which exact solutions 
are not known. 

In addition to providing a tool for investigating the dynamics of quantum spin 
models, our results also shed light on more general features regarding the role of closure 
conditions when truncating the BBGKY hierarchy: The particularly simple structure 
of spin- 1/2 lattice models facilitates analytic calculations beyond what can be achieved 
in continuum systems, and an understanding of the effect of approximation schemes 
is easier to attain. 

The article is structured as follows: In section[2]we introduce a fairly general class 
of long-range interacting quantum Heisenberg models with anisotropic interactions, 
subjected to an external magnetic field, and discuss its relation to the exactly 
solvable long-range Emch-Radin model. A short introduction to the quantum BBGKY 
hierarchy is given in section [3] with special emphasis paid to the role of closure 
conditions when truncating the hierarchy. In section l4~Tl we introduce a representation 
of the hierarchy, i.e. a certain choice of basis in the underlying Hilbert space, different 
from the conventional cluster expansion. In section 14. 3[ the BBGKY hierarchy is 
expressed in terms of this representation, and we find that it has the form of a second- 
order recursion relation. This structure turns out to be crucial for the derivation of 
our results. First, as discussed in section[5l the recursion relation suggests a scaling of 
the expansion coefficients and the corresponding time scales in powers of A 1 / 2 with 
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the system size N, implying a separation of time scales in the large system limit. 
The largest time scale is found to diverge proportionally to N 1 / 2 in the large system 
limit. As a consequence, equilibration is expected on a time scale that diverges in the 
thermodynamic limit, in agreement with what is known about the long-range Emch- 
Radin model. In section|51 we consider the long-range anisotropic quantum Heisenberg 
model for a special set of parameter values and initial conditions for which calculations 
are easier to perform. Under these restrictions, we can show analytically in section RTTl 
that closing the BBGKY hierarchy by neglecting £-spin correlations does never lead 
to equilibration, but gives rise to quasi-periodic time evolution with at most £/2 
independent frequencies. Moreover, in section f6. 21 we solve the untruncated recursion 
relation (full BBGKY hierarchy) in the large-A*" limit, observing a superexponential 
decay to equilibrium in rescaled time r = tN~ x / 2 . This behavior and the observed 
time scale are exact analytic results which are in a perfect agreement with the Emch- 
Radin model. A more detailed discussion of these results and their implications on 
the general considerations of section [3] is given in section [JJ and we summarize our 
findings in section [H 



2. Long-range anisotropic quantum spin model 

Two quantum spin models are introduced in this section. The first one, an anisotropic 
quantum Heisenberg model with Curie- Weiss type interactions, is the one actually 
studied in this work. The second one, named after Emch and Radin, is an Ising type 
quantum spin model in a longitudinal magnetic field, with two-body interactions that 
decay with the distance as a power law. For special choices of the parameters, both 
models agree, and the known exact results on the dynamics of the Emch-Radin model 
can be used to test the validity of our results regarding relaxation to equilibrium in 
the anisotropic quantum Heisenberg model. 



2.1. Curie-Weiss anisotropic quantum Heisenberg model 

Consider N identical spin-1/2 particles, attached to the sites of a d-dimensional lattice. 
The corresponding quantum dynamics takes place on the Hilbert space 

JV 

■*fr = (g)C?, (1) 
i=l 

where the C? are identical replicas of the two-dimensional Hilbert space of a single 
spin-1/2 particle. The unitary time evolution on M'n is generated by the iV-body 
Hamiltonian 

N N 

v »< ' X l " ( 2 ) 

i<j 

consisting of an on-site potential and a spin-spin interaction potential, 

H l = -J2 h a <, Vij =~ £ J ab <^, (3) 

ael a,b£l 

where erf denotes the Pauli operators at a lattice site i and 

l={x,y,z] (4) 
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is the set of component indices. The spin-spin interaction Vjj is of Curie- Weiss type, 
coupling each spin to every other on the lattice, and the strength of the coupling 
is determined by a 3 x 3 matrix J of coupling constants with matrix elements J ab . 
The 1/N prefactor in ([3]) is introduced to render the energy per spin finite in the 
thermodynamic limit. The long-range interactions, and the TV-dependent prefactor 
they are necessitating, can lead to peculiar properties. In equilibrium, it may happen 
that different statistical ensembles (like the microcanonical and the canonical one) are 
nonequivalent in the thermodynamic limit, and issues of this kind have been discussed 
in [19, 20J. In the present article, we will focus on peculiarities of the non-equilibrium 
behavior of this model. 

2.2. Emch-Radin model 

The Emch-Radin model [5J |5] has a Hamiltonian similar to ^ , but with a two-body 
potential 

V* = -pff^J (5) 

So in contrast to the anisotropic Heisenberg model's two-body potential ([3]), J zz is the 
only nonzero element of the coupling matrix J. Moreover, the spins are associated to 
the sites of a d-dimensional lattice, and the coupling between spins of and oj depends 
algebraically on their distance D(i, j) on the lattice, i.e. J zz oc D(i, j)~ a with some 
nonnegative exponent a. The on-site potential potential is of the same form as in ([3]), 
but with an external magnetic field h = (0,0, /i 2 ) pointing in z-direction. It follows 
from these definitions that the Hamilton operators of the anisotropic Heisenberg model 
and the Emch-Radin model have a special case in common: For a = the distance 
dependence of the Emch-Radin coupling J zz is eliminated and the Hamiltonian is 
identical to that of the Curie- Weiss anisotropic quantum Heisenberg model with 
J = diag(0, 0, J zz ) and h = (0, 0, h z ). 

The permissible initial states in the Emch-Radin model are restricted to density 
operators which are diagonal in the a x tensor product eigenbasis of J%n. Under this 
condition, an analytic expression can be obtained for the time evolution of expectation 
values of operators of the form 

N 

A = X>o? (6) 

i=l 

with real coefficients tjj [5]. In the case of short-range interactions, i.e. for exponents 
a > d, the time evolution of the expectation values (A) (t) was analyzed in the 
thermodynamic limit N — ¥ oo in [5] . The system was found to relax to equilibrium, and 
the process of relaxation was shown to be non-Markovian. The long-range scenario 
with $S a $S d, analyzed in [101 111) , was shown to display a remarkably different 
behavior: The time scale at which equilibration appears to occur was found to increase 
proportionally to N r with r = min{ 1/2,1 — a}. This finding implies a diverging 
equilibration time scale in the thermodynamic limit (see figure Q] for an illustration). 

The analytic results in [8] and |10[ [TT] are not easily extended to other parameter 
values, different initial conditions, or more general observables. We believe, however, 
that similar quasi-stationary behavior (i.e. long-lived non-equilibrium behavior persist- 
ing on a time scale that diverges with increasing system size) shows up under much 
more general conditions. The aim of this work is to employ quantum kinetic theory 
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Figure 1. Time evolution of the expectation value of the observable A for 
magnetic field h = 0, coupling matrix J = diag(0, 0, 1), and exponent a = 0. 
The plot was obtained by evaluating with Mathematica the analytic formula (9) 
of reference |10| for various system sizes N . The result is valid for arbitrary initial 
states and arbitrary coefficients in the expansion (J6j of A. The expectation 
value shows an apparent decay towards its zero equilibrium value, but on a time 
scale that depends strongly on the system size N (note the logarithmic time scale). 
The decay is only apparent as, on time scales much longer than shown, recurrent 
behavior (Loschmidt echos) occur due to the system's finite size. Similar behavior 
is observed for other values of a between zero and the lattice dimension d. 



to study the long-time dynamics of long-range interacting quantum spin systems of 
the Heisenberg type <j2j> and ((3]). For the moment, we restrict the analysis to Curie- 
Weiss-type potentials as in ((3]), but we expect qualitatively similar results to hold for 
algebraically decaying long-range potentials with exponents a between zero and the 
lattice dimension d. 



3. BBGKY hierarchy 



The density operator p^ of an iV-spin system is a self-adjoint, positive, trace-class 
operator, acting on the Hilbert space J^n- We use the normalization convention 

Tri...jvjOjv = 1, (7) 

where Tri. jy denotes a trace over all N factors of the tensor product Hilbert space 
((T|). The expectation value of an operator A with respect to pn is given by 

(A) = Ttl.jv ( Pn A) . (8) 

The time evolution of the density operator is governed by the von Neumann equation 

ihdtPN = [Hi...n, Pn] ■ (9) 

Reduced ^-particle density operators are derived from p^ by tracing out (JV — t) of 
the factors of J#at, 

Fi..j = Tt£+i...nPn, (10) 

where Tr^_|_i...jv denotes such a partial trace. The F\,..g are again density operators, 
i.e. self-adjoint, positive operators of trace-class. Since they are all derived from the 
same p^r, the reduced density operators are not independent, but satisfy a collection 
of consistency conditions of the form 

TnF 1 ... i = F 1 ... l -i (11) 
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(with the convention of F = 1) which we will refer to as the trace property of 

We will assume in the following that all reduced density operators Fx £ are 
invariant under ^-permutations of their indices from the set {1, . . . , N}. For £ = 1, for 
example, this means simply that 

F 1 =F 2 = ---=F N , (12) 

and for i = 2 this property amounts to 

F ij = F kl \/i,j,k,l€{l,...,N}, ijtj, k^l. (13) 

In the context of continuum (off-lattice) systems, permutation invariance is usually 
justified by the indistinguishability of particles. For the lattice systems investigated 
in the present article, the spin degrees of freedom can of course be distinguished 
by the lattice site they are attached to. Instead, in that context the assumption of 
permutation invariance amounts to assuming that the system is in a homogeneous 
state, i.e. different subsystems behave similarly. This assumption is expected to be 
justified in the Heisenberg model with ferromagnetic coupling J ab ^ 0, but will 
be violated in the presence of anti-ferromagnetism. Since the Hamiltonian @ and @ 
is permutation invariant as well, the homogeneity of an initial state will be preserved 
under time evolution. 

Under the assumption of permutation invariance, the Bogoliubov-Born-Green- 
Kirkwood-Yvon (BBGKY) hierarchy is obtained by applying the partial trace Tr^ +1 jv 
to the von Neumann equation ((§]). For a Hamiltonian of the form the hierarchy 
can be written as 

£ e e 

ihd t F 1 ... i = J2[H i ,F 1 ... e }+ \Vij,Fu.j] + (N-i)T[i + iY t \y^ + i 1 F 1 ..u+i]. (14) 

i— 1 i,J — 1 i=l 

i<j 

For a given number N of spins, this is a finite set of equations, and it is fully equivalent 
to the von Neumann equation ^ from which it was derived: A solution of the BBGKY 
hierarchy is equivalent to an exact solution for the density operator p^, from which 
the reduced density operators F\_ j can be obtained via equation (|f 01) . 

3.1. Closure conditions 

In the form (|14[) . the hierarchy is not yet particularly useful: Solving the full hierarchy 
is as difficult as solving the von Neumann equation (and therefore impossible for most 
cases of interest). For many physical problems, it turns out that n-particle correlations 
are naturally ordered by decreasing relevance, and only those with n ^ £ have to be 
considered (where £ is some small number, usually not larger than 4, depending on the 
system under investigation and the quantity of interest). It is therefore often sufficient 
to consider only the time evolution of the first £ reduced density operators. Hence a 
useful computational tool may be derived by truncating the hierarchy at the level of 
the £-spin reduced density operator. However, the achieved simplification of such a 
truncation comes at the expense of an approximation. 

Simply truncating the set (fT4"|) after the first £ equations results in an ill-defined 
problem: The £th equation, which determines the time evolution of iq...£, requires 
Pi...£fi as an input, but the equation determining iq...f+i has been eliminated by 
truncation. To obtain a well-defined system of £ equations, a closure condition has 
to be postulated between iq. ^+i and the lower order reduced density operators. We 
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refer to such a relation among i J i.„^+x and all of the {Fi...j} j—\ as an £th order closure 
condition. 

A frequently used truncation scheme of the BBGKY hierarchy is based on the so- 
called cluster expansion |21) . where correlation operators G\,,,i are defined according 
to the following scheme, 

F 12 = F l F 2 + G l2 , (15a) 

-F123 = FiF 2 F 3 + F1G23 + F2G31 + F3G12 + G123, (15o) 

-F1234 = F1F2F3F4 + F1F2G34 + F1F3G24 + F1F4G23 + F2F3G14 

+ F2F4G13 + F3F4G12 + G12G34 + G13G24 + G14G23 

+ -F1G234 + -F2G134 + F3G124 + -F4G123 + G1234, (15c) 

and so on. Based on this expansion, a straightforward way to close the BBGKY 
hierarchy at ^th order is to express (|14l) in terms of the correlation operators and set 

G 1 ... i +k = yi^k^N-e. (16) 

We shall refer to this approximation as the £th order correlation closure. 

Needless to say that the accuracy of the approximation will depend on the 
"quality" of the closure condition. Regardless of the details of the closure, it seems 
plausible to expect an improvement in the accuracy with increasing order I of the 
truncation. We have tested this expectation by numerically investigating the BBGKY 
hierarchy truncated by correlation closures of various orders. The density operators 
were expanded in the basis of Pauli matrices and the evolution was performed for the 
coefficients of the expansion (discussed in detail in section 0|) . The one-spin reduced 
density operator in this expansion reads 

^i = ^(li+EW), (17) 

where li denotes the identity operator on the one-spin Hilbert space. The real 
expansion coefficients /" are related to the mean spin expectation value, as 

N \ , N 



\ 1=1 / i=i 

The time evolution of the modulus 



l/ih^ + t/nH^) 2 (19) 

is displayed in figure [2] for correlation closures of orders I = 2, 3, 4. To summarize the 
simulations, the effect of increasing the order of the truncation is rather disastrous and 
the long-time evolution is badly predicted at either order: According to definition (fT0|) . 
density operators F\..ji have to be positive operators, and this property is conserved 
by the von Neumann equation ©. In the Pauli representation (|17p . positivity of F\ 
amounts to the condition 

l/i K 1 (20) 

which, as is evident from figure^ is violated for the correlation closures of order t 3 
after relatively short times. So the higher-order correlation closures not only fail to 
improve on the lower-order ones regarding the relaxation to equilibrium, but they 
even fail to preserve the basic features of the density operators. We will see in the 
forthcoming sections that the latter is an artefact of numerics, caused by the presence 
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Figure 2. Time evolution of the modulus of the coefficient vector f\ 

parameterizing the one-spin reduced density operator F%. The different curves 
correspond to correlation closures of orders i = 2, 3, and 4. The parameter values 
in the Hamiltonian © are h = (0,0, 1), J = diag(0, 0.04, 1), and N = 400, the 
initial conditions are chosen as /i(0) = (0.95,0,0) and gi(0) = for £ ^ 2 [see 
( 128 al l— (128 ctt l . The constant line at 1 indicates the boundary value of condition 
{20J at which F± ceases to be a positive operator. Remarkably, the larger the order 
I of the correlation closure, the earlier this property is violated. Using, instead of 
the original equations, averaged (Avg) time evolution equations as introduced in 
section [5] where the fast oscillations have been eliminated, the numerics improves 
significantly in the sense that F\ remains a positive operator for the times shown. 
The Gaussian plotted as a solid line indicates the behavior expected from the 
exact time evolution in the N — > oo limit. 

of oscillatory degrees of freedom with several fundamentally different time scales. This 
problem will be dealt with in section [5j where an averaging procedure is defined that 
eliminates the fast oscillations while correctly reproducing the evolution on the slow 
time scale. 

It is the main objective of the present paper to investigate and better understand 
the effect of the correlation closures at various orders, and in particular the question 
of whether and how this type of closure condition can correctly describe the relaxation 
to equilibrium in a long-range quantum spin system in the thermodynamic limit. 

4. Representation of the BBGKY hierarchy 

The derivation of the results reported in this article depends crucially on a particular 
expansion of the reduced ^-spin density operators in the basis of Pauli operators. 
Expressed in terms of this expansion, the BBGKY hierarchy reveals a recursive 
structure which is at the basis of the analytical results obtained. 

4-1- Definition of the expansion 

Inspecting the BBGKY hierarchy (fl4"|) . one might at a first glance get the impression 
that the reduced density operators iq...^ are coupled by a first order recursion, in 
the sense that the time evolution equation of F\_ j contains only one further reduced 
density operator, F\... However, this is not really true, as in fact the various Fi„j 
are dependent on each other through trace properties (fTTj) . The crucial idea behind 



Gaussian 

2" d order " " 

3 . order ~" 

4 order — 

3 order Avg 

4 order Avg 
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the expansion introduced in this section is to choose the expansion coefficients such 
that the trace property is automatically satisfied. Then the resulting coefficients are 
independent variables, on the basis of which the "true" structure of the equations can 
be studied. 

Expanding the first few reduced density operators in the Pauli basis 
{1, a x , o~ y , cr z }, one can verify by inspection that the following choice of expansion 
coefficients satisfies the trace properties (fTT|) . 



Fi =^(ii+E (21a) 

Fi a = - (l 12 + £ /i K + <tf) + E /2V^) , (21b) 



a&T a,b£T 



i \ A £ ab ( a b , a. b , _a b \ , \ ^ rabc a b c \ /oi \ 
+ ]Lj J 2 (< 7 lC r 2+ cr 2 0-3 + cr 3 cr l)+ 2^ ^ 0-^2^. (21c) 

a,b£T a,b,c£T 

We introduce several definitions to facilitate the general formulation. For the purpose 
of counting the particle (subscript) indices, we define, for 1 ^ n ^ I ^ TV, the set ty n {£) 
consisting of all n-element permutations (pi, . . . ,p n ) of particle labels pi £ {1, ■■■,£} 
such that pi < pj for all 1 ^ i < j ^ n (i.e. all sequences are strictly increasing). 
Considering for example n = 2, we have ^(-Q = {(1, 2), (1, 3), ... , (1, £), (2, 3), . . . , (£— 
!,£)}■ Moreover, we use the notation 



n<; ( 22 ) 

1=1 

with p £ ^Pn(^) and the component multi-index a = (oi, «2> . . . , a n ) with G I. 
With these definitions, the general pattern behind (j21a[) - (|21c[) can be captured by the 
formula 

l 

F 1 ... e = 2- £ J2 E fn E < (23) 

with the convention f§ = 1. 

By inspection of the last terms in (|21bl) and (|21c[) , it follows that f% b and f§ bc have 
to be symmetric with respect to permutations of their superscript indices in order to 
make sure that the resulting expressions for F\2 and -F123 are symmetric under particle 
exchange. As a result, not all components ff are independent variables. In particular, 

any element fi is equal to the element f% , where 

aim) = (x, ,x,y, ...,y,z,...,z) (24) 
m x elements 

is a permutation of a' with the components x, y, z arranged in contiguous blocks of 
m x , Tn y , and m z labels. Hence, all independent components of fi can be labelled 
uniquely by a triple of nonnegative integers 

m = {m x ,m yi m z ) 1 m, 5^ 0, (25) 

where m x + m y + m z = I. In the following, notation in terms of a and m = 771(a) will 
be considered as equivalent, ff = /"'. 
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The number of independent components of /J™ is equal to d(£) = (£+ l)(£ + 2)/2. 
Therefore the TV-spin density operator is described by 

N 1 

D(N) = E = q(N + 1)(N + 2)(JV + 3) - 1 (26) 
e=i 

independent real variables, a number substantially smaller than the 4 N real variables 
necessary to parameterize an arbitrary Hermitian operator on the TV-spin Hilbert space 
Ji? N given in ((TJ). 

4-2. Relation to the cluster expansion 

In order to implement the correlation closures (fl"6|). we need to express the correlation 
operators Gi,,,i in terms of the expansion coefficients /° 1 "' af . Substitution of the 
expansion (j2"3")) of F\...t into the cluster expansion (|15a|) - (|15c|) leads to 

G L ,=2-^ s ^ W)i (27) 

where gf are expressed in terms of ff by relations of the type 

If =9? + Kfx, (28a) 

fabc „aoc , fa rbc , fb fca , fc fab o fa fb rc /no/, 1 ! 

h =93 + /1/2 + /1/2 + /1/2 —*hhJi> \ 2 ° b ) 

fabcd „abcd i -fa rocd i x6 facd i t c fabd i rafoc 

/4 — 34 + /l /3 + / 1 / 3 + / 1 / 3 + J I J 3 

+ it ft + / 2 °7| d + ^"/^ 

- 2A b /i d /2 ac - 2/f/f/a* + e/f/f/f/i ■ (28c) 
The correlation closures of orders £ = 2,3,4 discussed in section [3~T1 are obtained by 
setting gi + i = 0. 

^.3. BBGKY in terms of the expansion coefficients 

Substituting the expansions (|2"3")) into the BBGKY hierarchy (fT4"|) , the time evolution 
equations for the coefficients fe can be obtained. After a cumbersome but rather 
straightforward calculation which is reported in |Appendix A| one finds the BBGKY 
hierarchy in terms of the coefficients fi to be 

\dtft = vm(ft) + H-(ft-i) + (1 - tyvf+Ut+i), (29) 
where the definition A = 1/N has been introduced. Moreover, we defined 



vM) = E h b J2^ bC fr a * +C > (30«) 
b.cei i=i 

£ 

«?-(/<-i)= "E E e* bc J ta '/r i 04+C "°'. (306) 



v a l+ {fi + i)=- E J bd E £ai&C /" + T +C+rf » (30c) 

bx,del 8=1 
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where e abc is the Levi-Civita symbol defined according to the convention e xyz = 1. 
With an abominable abuse of notation, 

a - a. t + c = (oi, . . . , aj_i,c, a i+ i, . . . ,a e ) (31a) 

in (|30a|) denotes the multi-index which is obtained from a by replacing its ith element 
by c. Similarly, 

a - a t + d - a 3 ■ — (ai, . . . ,Oi_i,d, Oi+i, . . . , aj+i, . . . , a*) (316) 

in (|306j) is derived from a by replacing the ith element by d and then deleting the jth 
entry of a, and 

a- a t +d + c= (oi, . . . , Oi_i, d, a 4+1 , . . . ,a e ,c) (31c) 

in (|30cp is obtained from a by replacing the ith element by d and then appending c. 

All the terms in viq and i>^± are linear in the coefficients fa and fe.±\. As 
inherited from the Hamiltonian, vio contains only the local field components h a , while 
V£± contain the spin interaction matrix J. Equation (|29p shows that the ith order 
coefficients jf are coupled to coefficients of orders I ± 1, but not to coefficients /" 
with \n — i\ > 1. It is this second-order recursion structure that is at the basis of the 
results derived in the following. 



5. The thermodynamic Limit 

The thermodynamic limit is a delicate issue in general, especially in the presence 
of long-range interactions. Recent work, mostly on classical systems, has taught us 
that complications increase further when one is interested in long-time asymptotics 
(or long-time averages) of long-range interacting systems, and their thermodynamic 
limit. The issue is illustrated by the time evolution of the expectation value (A)(t) 
in figure [T] where, for the Emch-Radin model, it was observed that relaxation to 
equilibrium takes place on a time scale that diverges with the system size N as N r , 
with an exponent r = min{l/2,l — a} > 0. As a consequence, depending on the 
order in which the long-time average and the large-system limit are taken, different 
limit values are obtained: Performing first the long-time averag^j] and then the large- 
system limit yields the equilibrium value (A) = 0, whereas the reverse order of the 
limits results in (A) = (A)(0). More generally, instead of taking one limit after the 
other, one can consider any path to infinity in the (N,t) plane and take both limits 
simultaneously along this path. It will depend on the physical question of interest 
which limiting procedure and which path is the suitable one. In this section we 
discuss one such path which, as we shall see, is appropriate for studying the approach 
to thermodynamic equilibrium of the long-range interacting spin system and (J3|). 

As always when performing a thermodynamic limit, it is essential to identify 
suitably defined quantities for which this limit is well-defined and nontrivial. With 
this aim in mind, consider the scaling transformations 

t =hr\- r /2, (32a) 

// = f^ s£ , (326) 

J For finite system sizes, the long-time limit is not well defined, as the time evolution is periodic 
and recurrences (Loschmidt echos) occur. The period of these recurrences, however, is exceedingly 
long for reasonably large systems, and the excursions away from the ensemble average are short-lived. 
For these reasons, the time-average of the expectation value (A)(t) has a well-defined long-time limit 
which, for large N, is close to the ensemble average. When the large-system limit is taken first, these 
problems do not occur, the long-time limit exists and is equal to the long-time average. 



Equilibration in long-range quantum spin systems from a BBGKY perspective 12 



where the parameters r and s are still to be determined. Substitution into (|29[) and 
multiplication by A~ r ~ s ^ yields the BBGKY hierarchy in the scaled variables, 

drf'e = A->o(/;) + A 1 - r - s ^_(/;„ 1 ) + X s - r (l - £X)v e+ (f' e+1 ). (33) 

The first term on the right-hand side of this equation depends, via veo, on the magnetic 
field components h a (i.e. on the on-site potential), but not on the spin-spin coupling 
J. The effect of such a constant field on the dynamics is known to be a simple spin 
precession, unrelated to the relaxation to equilibrium our analysis is focussing on. 
Moreover, as we will see later, its frequency turns out to be divergent on the time 
scale of equilibration. 

It is therefore convenient to hide the presence of the first term on the right-hand 
side of (|33|) by writing the BBGKY hierarchy (|14p in some kind of interaction picture. 
This amounts to absorbing the unitary time evolution, caused by the one-spin terms 
Hi in the Hamiltonian in the definition of a new set of reduced density operators, 



Ui. 



exp 



S5> 



Fi...£exp 



t5> 



The BBGKY hierarchy in the interaction picture is then given by 



, Ui...e + i 



where the two-body potential now is explicitly time-dependent, 



Vij(t) = exp 







}(H i + Hj )_ 


Vij exp 



(34) 



(35) 



(36) 



To compute the explicit form of Vij, note that Pauli operators on different lattice 
sites commute, and therefore the exponentials of sums in (|36|) can be factorized into 
a product of exponentials. Define the matrix B(b) as 

B{b)<j = e- lah ae lah 

= ct cos |2&|+ 26(6 -cr) sin 2 |6| + (erx S)sin|26|, (37) 

with b = (b x ,by,b z ), a = (o*,(T*,o*), \b\ = [(b x ) 2 + (bv) 2 + (b 2 ) 2 } 1 / 2 , and b = b/\b\. 
Applying ([57]) to the exponential factors in yields 

a,bel 

in a form that is equivariant to V12 in ([3]), with a transformed, time-dependent coupling 
matrix 

J = B(ht/h) T JB(ht/h). (39) 

Time appears only in the trigonometric functions in B, therefore B(ht/h) is periodic 
with period T = irh/ \h\. The elements of B(b) can be represented in the form 

r _j_ r .xx r xy _j_ c z r xz _ r y\ 



B(b) 

where c° = cos |26|, c % 



: +c x 

+ c zz 
2 



(40) 



b l sin \2b\, and c y = 2b l V sin 2 \b\ with i,j G 2 



Equilibration in long-range quantum spin systems from a BBGKY perspective 13 



Comparing the original BBGKY hierarchy (fT4"|) to the one in the interaction 
picture (|35[) . we observe that the single-spin term — whose diverging frequency might 
be bothering us in the thermodynamic limit — can be eliminated from the BBGKY 
hierarchy by passing to the interaction picture, at the expense of a time-periodic 
interaction potential (f3"6| . Next we introduce, entirely analogous to equation (|2"3"|) , an 
expansion of Ui...t in terms of coefficients u^, 

cu=2-TE< E < ( 41 ) 

Writing the BBGKY hierarchy in the interaction picture ([35| in terms of these 
expansion coefficients, we obtain 

a T «i = A 1 - r -'^_(t t J_ 1I r) + A'- r (l-a)^+(ui +lJ r). (42) 

where U( = u' e X se . The ^±(u^ ±1 ,r) terms are obtained from the corresponding ve± 
terms in (|30o[) and (|30cp by replacing J with 

J{t) = B{\- 1 / 2 hT/2) T JB(\- 1 / 2 hT/2) (43) 

as defined in (|3T))) . Inspection of (|4"2")l suggests that the leading-order A-dependent 
factors can be scaled to unitj{§| by choosing 1 — r — s = and s — r = 0. The solution 
r = s = 1/2 determines the scaling exponents in Q32ap and (|32fcl) , 

* =hrX- 1 / 2 /2 1 (44 a) 

/< = /;^ /2 , (446) 

and the BBGKY hierarchy simplifies to 

d T u' e = ve-(u' e _ 1: T) + (1 -£\)vt + (u' t+1 ,T). (45) 

The frequency of the modulation of J(r) in ([35]) has a singular A — > limit, 
which is a hallmark of the fundamentally different time scales involved: The fixed 
time scale of the single-spin interaction, contained in the term vio(f) in (|29|) in the 
(t, A r )-plane, appears as a variable frequency oscillation with a singular A = limit in 
the parameter space (r, A). Since this frequency diverges on the r-time scale we are 
interested in, it is reasonable to eliminate the high-frequency oscillations by applying 
an averaging procedure to equation (|45[) . To this purpose, we apply on both sides of 
(|45|) the finite-time averaging operator 

1 



T 



T + T 

dr 1 , (46) 



where T = TryX/\h\ is the period of J(r). For small A, ug is slowly varying over a 
period T, and by making use of the definitions (|30b[) and (|30c[) we can write 

6, eel i ,j=l 

(47) 



(1-A) ^ J 6rf ^ e °* 6c ^r +C+<i (r), 

b.c.del 



§ This argument is based on the assumption that v t and are of order unity in A. It will be made 
explicit for the special case discussed in section [til 
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(48) 



* Jo 

Performing the limit A — > (which implies T — > 0), we obtain the averaged BBGKY 
equations 



where V£± is obtained from the corresponding v#± by replacing J with J. 

Comparing the original (non-averaged) equations (|33p with the averaged ones in 
(|49l) . the outcome of the averaging procedure can be summarized as follows: 

(i) The presence of a single-spin potential Hi in the Hamiltonian leads to a fast 
oscillating term. 

(ii) Going to the interaction picture and averaging-out the fast oscillations, the 
averaged BBGKY hierarchy (|49|) is formally identical to the original one ([33)) 
in the absence of a single-spin potential, but with the original coupling matrix J 
replaced by the averaged one J defined in (|4"5|) . 

(iii) For the investigation of the long-time dynamics in the thermodynamic limit it is 
therefore sufficient to study the Hamiltonian @ in the absence of a single-spin 
potential Hi, thereby avoiding the problem of fast oscillations in the rescaled 
BBGKY hierarchy (|3%)). 

(iv) Our scaling ansatz (|32a[) and (|325[) suggests to define the thermodynamic limit 
as 1/N = A —> at fixed rescaled time r, and this suggestion will be further 
investigated in section [6] This limit corresponds to a non-trivial path to infinity 
in the (t, N) plane. 

Taking for example the special case discussed in section |3~T1 J = diag( J x , J v , J z ) and 
h = (0, 0, h z ), it is straightforward to calculate the averaged coupling matrix 



The above results then assert that, on the r-time scale defined in (|44a[) , equilibration 
of a system with coupling J and zero magnetic field will look identical to that of a 
system with the original couplings and a non-zero field in z-direction. For the special 
case studied in section |6] where an exact solution (without averaging) is feasible, these 
considerations will be nicely confirmed and illustrated. 

As a consequence of the elimination of the fast oscillations, the averaged time 
evolution equations (|49p are much more suitable for numerical computations, yielding 
more stable results. Comparing the results with and without averaging, we observe 
that the reduced one-spin density operator retains positivity over the period of time 
shown in figure[2j Moreover, the solution appears periodic, which is in agreement with 
the analytic prediction for the special case investigated in section [6] 

6. Case study 

In this section we derive, for the choice of parameters 



d T u' e = vi-(u' t _i) +Vi + (u' e+1 ), 



(49) 



J = diag {{J x + J*)/2, (J x + J y )/2, J z ) . 



(50) 



J = diag(J , J ,J Z ), 



h = (0,0, h z ), 



(51) 
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the time-evolution of certain coefficients fg in the BBGKY hierarchy This is 

achieved for correlation closures of arbitrary order, and also for the full (untruncated) 
hierarchy. As for initial conditions, we set 

(/i,/i tf ,/f)(0) = («*>*">**). H<1. (52) 
The parameters (fSTj) and the initial conditions ([52]) are more general than the ones 
used in the Emch-Radin model (where J = and the initial density operator is 
diagonal in the a x eigenbasis, implying f\ = /-f = 0). For the higher order coefficients 
we require 

tf(0) = s? = W>1. (53) 

A study of more general initial conditions and parameters will be presented elsewhere. 
The results we will obtain in this section are found to be in agreement with the 
known exact solution of the Emch-Radin model. This can be seen as an a posteriori 
justification of the scaling assumptions (j44a[) - (|44b[) and of the limiting procedure 
proposed in section [5] Apart from this reassuring result, the derivation is interesting 
to follow as it resorts to a number of beautiful mathematical concepts, including 
recurrence relations, continued fractions, and orthogonal polynomials. 



6.1. Finite system size 

Unless stated otherwise, we will drop from now on the prime in the rescaled coefficients 
(|326p . writing / instead of /'. We start by analyzing the hierarchy (|3"3"|) , and the terms 
vio and vi± contained therein, for the above parameter values and initial conditions 
(|5ip - (|53|) . In the following, whenever convenient, we will label the coefficients /J™ by 
a triple m = (m x ,m y ,m z ) as introduced in section f4. II Considering m = (0,0, £), we 
find veo = and 

d T f^ = t{J^ J^/itf'- 1 ' = 0, (54) 
from which we can conclude that 

/f°^( r ) = /f°>>) = s (W) (55) 

for all I. In particular, for I = 1 we have 

/f (r) = /? (0) = s\ (56) 



From (|55|) , we also note that the initial conditions (|53|) may be generalized to nonzero 
s (o,o,-0 a t no additional expense. Next we discuss the time evolution of coefficients /™ 
for two families of superscript indices, m = (1, 0,^ — 1) and m = (0, 1), 

d T f^ ^ = A-V^yf M-D _ {i _ 1)A7 (o ,14-2) _ (1 _ xe)Kf^/\ (57a) 
d T fp l ' l - 1] = \~V*h' + (£- l)Kf£°/-V + (1 _ X£)Kf^ e \ (57b) 

where K — J z — J , The regularity behind this set of differential equations can be 
captured with the help of a few definitions. Define the sequences of component labels 

a x = (x, yz, xzz, yzzz, xzzzz, . . .), (58a) 
a y = (y, xz, yzz, xzzz, yzzzz, . . .), (586) 

and denote by a^j and the nth element of the corresponding sequence, 1 ^ n ^ N. 
Furthermore, define a sequence of complex functions 

u„(r)=e iA " 1/2 ^[(-l)L"/ 2 J/5(r)+i(-l)L(«- 1 )/ 2 J/5(r)], (59) 
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where \n/2\ denotes the integral part of n/2. With these definitions, the BBGKY 
equations (|57ap and (|576p can be expressed in terms of the functions u„, 

d T u n = K(n - l)u„_i - K(l - n\)u n+1 , I ^n^N. (60) 

This equation may be further simplified by another scaling of time r' = t/K , yielding 

d T u n — (n - l)u„-i - (1 - n\)u n+ i, 1 n s$ TV, (61) 

where we have dropped the prime from r' for the sake of a simpler notation. 
Translating the initial conditions (|52p and (|53|) to our new variables, we obtain 

u(0) = (s x +is y ,0,...,0). (62) 

Inspecting equation (f5U|) . we observe that indeed the singular term X~ 1 ' 2 veo in ([55)1 
results in a high-frequency oscillation superimposed on the much slower dynamics 
induced by the spin-spin interactions. Incorporating this oscillation into the definition 
of u corresponds to a change of variables to the co-moving frame, i.e. to the interaction 
picture defined in section [5j Furthermore, (|59[) justifies the assumption, mentioned in 
a footnote in section [5J that r>e- and vi+ are of order unity in A. This gives support 
to the reasoning behind the choice of the exponents r and s in the scaling ansatz 
(132 al) (|32fe[) . However, it also implies that, on the r-time scale, no truncation of the 
BBGKY hierarchy can be justified by the smallness of the parameter A. 

Finally, the particular choices of parameter values (|5ip and initial conditions (|52p 
and ([53)) led to a simplification of the BBGKY hierarchy. The simplification becomes 
even more pronounced if one focusses, as we did in the preceding paragraph, onto a 
certain subset of the coefficients ff in the expansion (|23|) . For example, in order to 
compute expectation values of single-particle observables, it is sufficient to determine 
the time-evolution of ff , ff, and /f . From (|6"Tj) , we can read off that /f + iff is 
determined as a solution of TV equations with TV independent variables, which is a 
tiny fraction of the total number of degrees of freedom D(N) of the full BBGKY 
hierarchy [D(N) ~ TV 3 /6 as given by (|26p j. It is important to note that this reduction 
of complexity is obtained without truncation of the hierarchy, but is a consequence of 
the decoupling of certain subsets of the expansion coefficients introduced in (|2"5)) . 

The initial value problem specified by (|6"Tj) and (|f)2"j) can be solved in Laplace space. 
The Laplace transform of a function u of a real variable r is a complex function u of 
a complex variable z. defined by 

/•OO 

u(z) = JSf [/] = / u(T) C - ZT dT. (63) 
Jo 

The inverse of the Laplace transform is defined as 

u{t) = ££~\u\ = — [ u(z)c ZT dz, (64) 
2tti J 1 

where 7 denotes the so-called Bromwich contour in the complex z plane, which runs 
from — ioo to +ioo and stays to the right of all of the poles of u{z). The Laplace 
transformation of a derivative is 

if[d r u(r)] = zu{z) - u(0). (65) 

By virtue of this rule, we can write (|6ip in Laplace space, 

zu n (z) = u„(0) + (n - l)u„_i(z) - (1 - nA)«„ + i(z). (66) 
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From equation (|66p with n = 1 we obtairfjj 

Ul (0)/ Ul (z) = z + (1 - X)[u 2 {z)/ Ul (z)}, (67) 
and for n > 1 it can be rearranged to 

z = (n - l)[u n -i(z)/u n {z) - (1 - n\)[u n+1 (z)/u n (z)}. 
We can then write the ratio of subsequent 

u n (z) n — l 



u n -i(z) z + (1 - n\)[u n+ i(z) I u n (z)] 
Applying this expression iteratively, we arrive at 



(69) 



with 



where 



«i(z) = (s +is v ) — — (70) 

z+ z+ z 



(3 n = n(l - nX), l^n^N, (71) 
ai a,2 a\ 



bi+b 2 + bl 



a-2 



(72) 



b 2 H 

is a standard notation for continued fractions. 

From Mi, the coefficients u n with n > 1 are obtained via the equation (|66p . For 
a finite number A = 1/A of spins, the solution (|70[) contains a finite number of terms, 
terminating at (3^ = N(l — NX) = 0. Hence we can write 

u 1 ( z ) = { s*+is«)^fl (73) 
B N (z) 

where A^(z) and Bm{z) are polynomials in z of degree TV — 1 and N, respectively. 

From ui we reconstruct the coefficients ff and ff in the time domain via the 
inverse Laplace transform (fM)) . 

/i T (r)+i/ 1 S (r)=e^ 1/2 ^^- 1 [^i(z)](r). (74) 

Properties of ([74]) depend crucially on the properties of Jz? -1 [ui(z)], and therefore on 
the structure of the singularities of u\(z). Owing to the form of the expression ([75)1 . 
these singularities are poles, originating from the zeros of Bm(z). We shall study these 
zeros in the following and discuss the effect of correlation closures of various orders I 
on the solutions. For this purpose, consider what is called the ^th convergent of the 
continued fraction ([70)1 . 



u?\z) = (s* + i S y)l- A...^zi (75) 



z+ z+ z 



defined as the truncation of (fT0"|) after € terms. The same u\ could also be obtained as 
a solution of a truncated version of the recurrence relation (|66[) where iie+k = for all 
1 ^ k ^ N — £. For this reason, the convergent ([75| is closely related to the £th order 
correlation closure of the BBGKY hierarchy, as given by (|28ap - ()28cp with gt+k = 0. 
This can be seen by applying the correlation closure condition to the coefficients 
fill with component indices from the sets (|58a[) or (|58 6[) : Because of the (n — l)-fold 



[| We disregard the problem of zero denominators in the following, as it can be circumvented and 
does not cause serious problems. 
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occurrence of z-indices in the elements and a^, the cluster expansions (|28a|) - (|28c|l 

a v 

are particularly simple. The coefficient / 4 4 for example can be written as 

f < = = n _ 6/f/r + 6(/f) 3j +3/r [/r _ 2{f z )2 ] + 3 /rvf+ffr « (?6) 

In third order correlation closure, i.e. setting g% zzz = 0, the right-hand side of ([75)) 
consists of only two kinds of coefficients f£: either with a € a^Ua", or with a = z . . . z. 
The latter ones, as we have shown in (|55[) . do not vary in time, and f% zzz is therefore 
linear in the time-dependent coefficients f£ v . Translating this observation via 
definition (|59|) into the interaction picture, we find the closure condition 

3 

ua{t) = ^c„u„(r), (77) 

n=l 

where the c„ are determined by the time-independent coefficients f™ with m = 
(0,0, k). The same kind of pattern emerges also for correlation closures of arbitrary 
order £, where 

l 

u t+ i(r) = ^ CnU n {r) (78) 

n=l 

is again linear in the time-dependent coefficients u n (r). In the special case of initial 
conditions s(°'°> fc ) = /(°>°>' c ) = 0, we have c„ = and the Hih order correlation closure 
is equivalent to setting u„ = for n > I. This condition translates into the condition 
ui+k = for the coefficients in Laplace space, and ([75)) is therefore the solution of 
the ith order correlation closure for initial conditions given by ([55| . In case of more 
general initial conditions with s(°'°' fc ) ^ 0, the correlation closure amounts to a linear 
relation between it^+i and all the lower order coefficients. 

Having understood the relation between correlation closures (|16[) and convergents 
(|75p. we proceed with the analysis of the latter. Similar to the full continued fraction 
in ()73[) . also the convergent can be written as a ratio of polynomials Ag and Bi of 
degrees I — 1, respectively £, 

= & + (79) 

From ([75)1 . Wallis-type [55] second-order recursion relations for B n and A n can be 
derived, 

zA n (z) = PnAn-^z) - A n+1 (z), (80a) 
zB n (z) = p n B n _ x {z) - B n+1 (z), (80b) 

where j3 n is defined as in ()7ip . Although of the same form, the two recursion relations 
(|80ap and (|80 &|) have different initial conditions, _B_i = 0, Bq = 1, and Aq = 0, 
A\ = — 1. Anticipating that the zeros of i? n and A„ are imaginary, we make a 
coordinate transformation x = — \z that brings A„ and B n into the conventional real 
form. In these new variables, the convergent can be written as 

flW (ia;)= _ i(8 - + ia , ) ?yM > (81) 

'(x) 

with polynomials 

pW(x) = -i n A n+1 (ix), pW(s) = i n fl„(Lr). (82) 
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From (|80ap and (|806[) . recursion relations for the new polynomials are obtained, 

xp^(x)=^ n+j p^ ) _ 1 (x)+p^ +1 (x), j =0,1, (83) 

with initial conditions p^\ = and Pq =1. It follows from Favard's theorem [55] 
that the sequences of polynomials {pn }„=o generated by (|83|) are positive definite 
and orthogonal, in the sense that there exists a nonnegative weight function w$ (x) 
such that, for any 0^n<m^N — j, 

plP(x)piiHx)w^(x)dx = 0. (84) 

It follows from the positive definiteness and orthogonality of the pffl that all £ zeros x^ k 
ofpg (x) are real, simple, and isolated. Most importantly, the following decomposition 
formula holds, 

(o), ; -2^ a ._ a . / ' aek = ^w ( — ? (8o) 
p e {x) fc =i x lk p'yixik) 

where the prime denotes a derivative. 

Equation (|85p can be used to perform the integral in the inverse Laplace transform 
(|64p . Since we have to undo the coordinate transform z = ix, the Bromwich integration 
contour 7 in the complex z-plane has to be modified to 7' in the complex x-plane, 
running parallel to the real axis and below all the zeros of pf'(x), i.e. below the real 
axis. The exponential e 1XT in introduces a damping factor if 5x > 0, and the 
integration path 7' can therefore be closed in the upper half plane by a path 7' without 
changing the value of the integral. As a result, the path encloses all £ zeros of pf (x) 
and we can apply the residue theorem, 

M, x _ -i(s x + i S y) f P ( A( x ) e 



(t) = 2^ / , — WT~ dx 

Z7r J 7 P e (x) 

= («■ + E r*. , (4^tt) eixe3T (so) 

= (s x +i S y )^2a ej e ix ^ T . 

i=i 

By inspection of (|83p , it is not too difficult to observe that the polynomials generated 
by this recursion relation are alternatingly odd and even functions. Their zeros are 
therefore distributed symmetrically around the origin. Sorting the zeros in increasing 
order, (|86p can be simplified to 

e/2 

u[ e \ T ) = 2{s x + is y ) a H cos (^-t), (87) 
3=1 

where we have assumed that £ is even. (If I is odd, a constant term must be added to 
the sum.) 

Equation (|87p is the main result of this section. Since it consists of a finite sum 
of oscillatory terms, the solution u\ (t) is either periodic or quasi-periodic in time. 
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Figure 3. For parameter values h = 0, J = diag(0, 0, 1), and N = 1/A = 10, the 

real part 3Z(u\ ) = /f is plotted as a function of rescaled time r for various orders 
I = 2,3, 4, 5, 9 of the correlation closure. Recurrent behavior is evident from the 
plot, both for the full solution and for the truncations of various orders. 



For example, by setting A = in the polynomials p/_j with I = 2, 3, 4, we hnd 

,( 2 ) 



u\^(t) = (s x + is v )cos: 



(88 a) 



«w (r) 



{s x + \s y ) 
3 

{s x + i s y) 



— 2 | i 



(V6 + 2) 



V61 



io) 



The behavior of the solutions in (|57|) is illustrated in figure [3] for various orders £ of 
the truncation of the continued fraction (|75|) . Initial conditions satisfying s z = were 
chosen, and therefore the truncations correspond to £th order correlation closures. 
The recurrences, or Loschmidt echos, of the solutions, expected from the result in 
(|87p . are evident in the plots, both for the full solution for N = 10 spins as well as for 
the truncations at orders i < N. 



6.2. Thermodynamic limit 

Having observed that, for any finite system size and/or finite order correlation closure, 
the dynamics is periodic or quasi-periodic, we next want to investigate the effect of 
the thermodynamic limit A — > 0. The reader is reminded that this limit, due to the A- 
scaling of time (|32ap , corresponds to simultaneously taking long-time and large-system 
limits along a nontrivial path in the (t, N) plane. For the special case studied in this 
section |6l this limit can be performed within the framework developed in section 16.11 
by studying p$_i/p$ m t ne limit A — > 0. Here, however, we prefer to take a shortcut 
and investigate the A = case directly in the time domain. Indeed, setting A = in 
(|6ip . we obtain 

u„+i(t,0) = (n - l)u n _i(r,0) - d T u n (r,0), (89) 
which is solved by 

u„(r, 0) = (s x + isvy 1 - 1 ^ 2 / 2 . (90) 
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Translating this result to the original coefficients f n n by ([55]) . we encounter a term 

■ , — 1/2 i z 

e r , oscillating at infinite frequency in the A — > limit. This oscillation is 

irrelevant for the relaxation to equilibrium, as it happens on a fundamentally different 
timescale. Therefore we will disregard this oscillating term, obtaining 

/? (r, 0) + i/f(r, 0) - («» + i^e"^ 2 (91) 

for n = 1. 

The solutions (f9TH) or (|9"Tj) display Gaussian (superexponential) relaxation as a 
function of rescaled time r, and no recurrent behavior is observed. This is in stark 
contrast to the periodic or quasi-periodic dynamics (|88a|) we found, as illustrated in 
figure |3l for finite system sizes. Moreover, the result in (pJl) is in exact agreement 
with the N — > oo limit of the Emch-Radin model with exponent a = as reported in 
equation (32) of [TT] , 

7. Discussion of the results 

The special case treated in section [5] illustrates and confirms some of the aspects 
that have been discussed in sections 13.11 and [S] in greater generality, and it also can 
substantiate some of assumptions made. In particular, we would like to discuss the 
following aspects. 

7.1. Symmetry-induced causal relations 

We have seen that, for the special case of section [HI not all of the expansion coefficients 
(or /") influence the time evolution of each other. Instead, we observed that the 
time evolution of, say, /f is affected by only a small subset of the other coefficients. 
This decoupling, which is a property of the structure of BBGKY hierarchy of equations, 
is particularly pronounced in the case discussed in section [5] but, as will be shown in 
a forthcoming paper |23) . similar (but more involved) relations hold true for general 
parameter values of the model. 

7.2. Shortcomings of the correlation closure 

For the special values of parameters and initial conditions of section [BJ we found that, 
for any I 1, the correlation closure of order I [as defined in (p~6])] gives rise to periodic 
or quasi-periodic behavior. There are, of course, many other types of closures of the 
BBGKY hierarchy one might consider, and the question is whether others might be 
more suitable for studying the relaxation to equilibrium we are interested in. 

One way to improve on the correlation closure might consist in taking into 
consideration the symmetries imposed by the BBGKY hierarchy as mentioned in 
section 17.11 Since symmetries are respected by the exact dynamics, one may hope 
that a closure which takes into account this structure might be superior to the one 
that does not. Preliminary numerical studies of ours seem to be in favor of this 
conjecture. However, it is clear from these results that for our main objective, i.e. 
the study of relaxation to equilibrium, such a modified closure does at best lead to 
marginal improvements and we will hence not pursue this aspect further in the present 
article. 

A kinetic theory appropriate for studying relaxation to equilibrium is expected 
to require a more complicated closure relation in the form of a collision integral. Such 
a theory will be developed in a forthcoming paper. 
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7.3. Thermodynamic limit and the definition of reseated time r 

A crucial step for obtaining a nontrivial long-time asymptotics was the definition of a 
suitable kind of thermodynamic limit as elaborated on in section O This limit reflects 
in the choice of the exponent r when defining the rescaled time variable r in (|32ap 
0. In the general considerations in section [3 the choice r = 1/2 was based on the 
scaling properties in A of the BBGKY hierarchy (|42|) in the coefficient expansion. With 
this choice, and for the special case of section El we were able to derive, in (|9Tj) , the 
exact Gaussian relaxation to equilibrium, which gives support not only to the choice 
r = 1/2, but also to the guiding principle based on the scaling properties of (|42D . 

8. Conclusions 

We have studied, via a BBGKY-type approach, the time evolution of ^-spin reduced 
density operators for Heisenberg-type quantum spin models with Curie Weiss- 

type long-range interactions. Our analysis is based on a particular expansion of F\...i 
in terms of coefficients ff, as introduced in (|23p. which casts the BBGKY hierarchy 
into the form of a second-order recursion (|29l) . 

Originally our study was motivated by the observation of quasi-stationary 
behavior in the long-range Emch-Radin model, i.e. the fact that relaxation to 
equilibrium takes place on a time scale which diverges with the system size. As a 
consequence of this diverging time scale, for studying the long-time asymptotics of 
quantum spin models in the thermodynamic limit it is therefore necessary to consider 
a suitably defined thermodynamic limit where, with increasing system size TV, time 
is scaled appropriately such that nontrivial dynamics can be observed. Remarkably, 
the structure of the BBGKY hierarchy (|%2")l when expressed in terms of the expansion 
coefficients ff suggests a definition of rescaled time r = 2t\ r /h which, as turns out, 
leads to a nontrivial exact result and reproduces the Gaussian relaxation (|90|) of the 
Emch-Radin model. 

When dealing with the BBGKY hierarchy P2"| in rescaled time, we noticed that 
the on-site potential (or single-spin potential) Hi in the Hamiltonian @ gives rise to 
an oscillating term whose frequency diverges on the r-time scale in the thermodynamic 
limit. Since this time scale of oscillation is strongly (infinitely) separated from the time 
scale of relaxation to equilibrium, an averaging procedure is introduced to eliminate 
the high-frequency dynamics. The averaged BBGKY hierarchy leads to a well- 
defined A — > limit of the hierarchy, but also to a significant improvement of numerical 
results as shown in figure [21 

A general BBGKY hierarchy cannot be solved exactly, as this would correspond 
to an exact solution of the full TV-spin problem. The problem becomes more tractable 
by truncating the hierarchy. The resulting set of equations is ill-defined, but it 
can be turned into a well-defined one by postulating a closure condition. For the 
long-range spin models discussed in the present work, we have defined what we call 
the correlation closure of £th order in (fT6| . and the effect of these closures on the 
long-time dynamics has been discussed. For the special parameter values and initial 
conditions considered in section[5J we have shown analytically that closing the BBGKY 
hierarchy by neglecting £-spin correlations does never lead to equilibration, but gives 
rise to quasi-periodic time evolution with at most £/2 independent frequencies, as is 

If Without rescaling of time, the time evolution of the variables u e in the interaction picture would 
just be constant in the thermodynamic limit, with no sign of relaxation to equilibrium 1111 . 
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evident from (|87|) . We must therefore conclude that, in order to construct a kinetic 
theory appropriate for studying relaxation to equilibrium, a more complicated closure 
relation, presumably in the form of a collision integral, will be needed. Once such 
a theory is properly benchmarked against the exact results available for the Emch- 
Radin model, it should allow us to study the approach to equilibrium in non-integrable 
generalizations of the Emch-Radin model, or non-integrable long-range variants of the 
Heisenberg model for which exact solutions do not exist. 

In addition to providing a tool for investigating the dynamics of quantum spin 
models, our results also shed light on more general features regarding the role of closure 
conditions when truncating the BBGKY hierarchy: The particularly simple structure 
of spin-1/2 lattice models facilitates analytic calculations beyond what can be achieved 
in continuum systems, and an understanding of the effect of approximation schemes 
is easier to attain. We tend to believe that, at least on a qualitative level, several of 
our observations apply not only to spin-1/2 systems, but should hold more generally 
for closed quantum systems on finite-dimensional Hilbert spaces. 

We want to conclude with a comment on an interesting related work by Sciolla 
and Biroli [24] which also deals with the dynamics of Curie- Weiss-type (or completely 
connected) quantum systems, i.e. with Hamiltonians of type ([2|), but with arbitrary 
Hi and Vij. Their analysis yields exact analytical results in the thermodynamic 
limit for a much larger class of Hamiltonians, but is more restrictive with respect 
to initial conditions, allowing only wave packets whose width shrinks to zero in the 
thermodynamic limit. Moreover, their analysis is bound to fail on the rescaled time 
scales t studied in the present work, and is therefore unsuitable for investigating 
thermalization. In a future work we will study thermalization in rescaled time for 
more general models by approximative methods, and the results by Sciolla and Biroli 
should prove useful for benchmarking the short-time dynamical behavior. 



Appendix A. Derivation of equation (1291) 

In the course of the proof, we will refer to the following elementary identities for Pauli 
operators on lattice sites i = 1,2. 

[a?,af] =2i]T e Qb X, (A.l) 

C 

\„o,—U —b—V~\ r\- \ ^ (sruv-dbc—C , ra6_u^c^.c\ /a n\ 

[<r 1 (J 2 ,a l a 2 \ =2l Z^\ d £ a i+ d £ a 2) j \ A - 2 ) 

c 

Tr 2 [of erj , fal] = 4i<5™ ]T e abc ^, (A.3) 

c 

with component (superscript) indices a, b, c,u,v £ 1. We use the notation S ab for the 
Kronecker tensor, and e abc for the Levi-Civita tensor (with e xyz = 1). In the identities 
involving Pauli operators on different lattice sites, the identity operator is implied on 
each single-spin Hilbert space Hi not acted upon by any of the Pauli operators [as on 
the right-hand side of (|A.2|) ]. 

The starting point of the proof is the BBGKY hierarchy in the form (|14p with a 
fixed value of £. Inserting the expansion ([23|) into the left-hand side of (fT4|) , we obtain 
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where T = {x, y, x} denotes the set of component indices. Since the operators cr" are 
linearly independent among each other, we can separately equate their coefficients in 
(|A.4|) to those which result from expanding also the right-hand side of (|14[) in terms 
of (|23p . We will in the following set n — £, as this will be sufficient in order to obtain 
time evolution equations for all coefficients Since = {(!,...,£)} contains 

only a single element, the sum over p in (|A.4|) consists only of one term, er^ The 
rest of this appendix will therefore be devoted to computing, from an expansion of the 
right-hand side of (fT4|) . all the terms proportional to (7% e y and then equate them 
to 

2- l ihB t fta a {1 ^ ) (A.5) 

in order to obtain the time evolution equation of We will discuss the three terms 
on the right-hand side of (fl"4"|) separately, showing that 

(1X51) = (TXT3l) + (1X201) + (1X231) . (A.6) 

We start by inserting the expansion (f2Uf into the first term on the right-hand side 
of {HI, yielding 

j2[H i ,F 1 ...t] = -2-*Y,h b f: Y,fn E ( A - 7 ) 

i=i bex i.n=iaex n pesPnW 

For the commutator in (|A.7|) we find 

3=1 3=1 c£X 

where (|A. 1|) has been used. Here, p — and a — aj denote the sequences obtained 
from p and a by deleting their jth elements, and 

a — a,j + c = (oi, . . ., Oj-i, c^aj+i, . . . ,a„) (A. 9) 

jth element 

is the sequence where the jth element has been replaced by c. We observe that, for a 
given p = (pi, . . . ,p n ) with n elements, the commutator in (|A.8|) is again a product of 
n Pauli operators acting on different lattice sites. Since, as explained above, we can 
restrict our attention to the terms proportional to a fixed it is sufficient to 

consider the terms with n = £ in (|A.7|) . 

2^i E ^ E E E /fe a3&c <" aj+c J (A.10) 

where the expression (|A.8|) for the commutator has been inserted. Since tye{£) = 
{(!,...,£)} consists of only a single element, the sum over p disappears. Moreover, 
for the same reason, we have <5 !iPj = <5 ? J , which allows us to execute the sum over j, 
yielding 

£ 

2 i -h J2 h b E ft E £a * bc<T l wr ( A - n ) 

6,c£l aGl ? »=1 

Swapping names of the summation indices a,i and c and making use of the cyclic 
property of the Levi-Civita tensor, we can rewrite (|A.11|) as 

-^E^ E h b J2e a ' bc fr a ' +c - (A.12) 

aGl f 6,cGX t=l 
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Owing to the mutual independence of the cr^, we now can pick, for a given a, the term 

i 

- 2 1 - £ i<r^ /) J2 h b Y / e a ' bc f^ +c (A.13) 

b,c£X i=l 

which gives the first contribution to the time evolution equation (|A.6|) for the 
coefficient f%. 

To deal with the second term on the right-hand side of (|14[) , we again replace the 
i!-spin reduced density operator F\„j by the expansion (|23p , obtaining 

{: y«, = -E^EEE/: E (a.m) 

j<i=l b.cel i,j=ln=la6l» p6<p„(£) 

Under the condition i j, we can write 

n 

\ „b „c _al \ ^ r r -.a — a k — ai \bn a k an 

fcj=l 



(A.15) 



fc=i 



for the commutator in (|A.14[) , where we have used the short- hand notation Sj& p = 
rim=i — 3j,Pm,)- The first term on the right-hand side accounts for the cases where 
both, i and j, have a counterpart in cr^ , the second term for when only one of them has. 
(If neither i nor j have a counterpart in cr° then the commutator is zero.) Applying 
()A.2j) to the commutator in the first sum in (|A.15j) . all summands will consist of 
products of n — 1 Pauli operators. Since n ^ £, none of these terms can give a 
contribution to (|A.14j) which is proportional to a'K e y and the first sum in (jA~T5]) 
can therefore be neglected. Applying (jA.lj) to the commutator in the second sum in 
(|A.15[) . however, results in products of n + 1 Pauli operators. Hence, inserting (|A.15[) 
into (|A~T4| . we can restrict the summation to the n = £ — 1 terms, as only those may 
lead to contributions proportional to er^ ^ , 



^ jbc 

E wn E ft-i E E ; 



-i 

E E a p-pt 

b.cei aei'- 1 p&$ t -^{l) i-.3=i fe=l (A.16) 

The terms in the round brackets in (|A. 16|) are symmetric to each other under the 
exchange of i o j, and their contributions to the overall sum are therefore identical. 
Hence we can write 

E E ft-i E E Y,^^ akbd < + r +d+ ^ (a-17) 

b,c,del atl 1 - 1 petyt-lit) i,3=l k=l 

where (| A. 1[) has been applied to the commutators in (|A.16|) . The set ^_i(£) consists 
of all sequences where precisely one of the numbers 1, ... ,1 is omitted. For example, 
for I = 4, we have q3 3 (4) = {(1,2,3), (1,2,4), (1,3,4), (2,3,4)}. Due to the constraint 
Sji^p in (|A.17[) . j has to equal this omitted number for all non- vanishing contributions 
to the sum. We can therefore combine every p € ^_i(f) together with j £ p into a 
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new summation index p' £ y$e(£), and then use the sum over j to count through the 
previously omitted elements of the sequence. With these new summation indices, we 
can rewrite (|A.17[) in the form 

2 i-l. t 



"E E E e aM J ba *ft? E °7 ai+d ' (A.18) 



N 

b,d£l 1,3=1 ael* p'&#i(t) 

where the former component index c has now been included into the multi-index a. 
Observing that 9(it(£) = {(!,...,£)} consists of only a single element and after some 
renaming of indices, we can write (|A.18|) as 

" E E E e^ d J^ftr d ~ a3 - (A.19) 



N 

a£X e b,d£X i.j = l 



Owing to the mutual independence of the cr^, we now can pick, for a given a, the term 

J- £ a,bdjba ]f a-a l+ d-a^ 



--^r»(i,..^) 

b,d£l i,p=l 

which gives the second contribution to the time evolution equation (|A.6|I for the 
coefficient Similar to the notation introduced in (|A.9j) . the multi- index 

a — Oi + d- a,j = (ax, . . . ,Oi_i,d, Ot+i, . . . , aj-x, a j+ x, . ..,at) (A.21) 

in (|A.20[) is derived from a by replacing the ith element by d and then deleting the 
jt\x entry of a. 

The contribution from the third term on the right-hand side of (|14[) is derived in a 
very similar way, so we will only sketch the calculation here. Replacing the {I + l)-spin 
reduced density operator i<i..^+i by the expansion (|23|) . we obtain 

l 

{N-t)j2Tr e +i [V^ +1 ,F 1M ] 

1=1 i + i t (A. 22) 

= -^E j6c EE/n E E^W^i.^]- 

6,cez n=iagi" peq3„ (^+i) 4=i 

Applying ()A.3|) to the partial trace of the commutator in (|A.22[) one finds, similarly 
to the discussion of the second term, that only summands with n = I + 1 can lead 
to terms proportional to rr^ « in (|A.22[) . As a consequence, the sum over p again 
consists of only a single term and, after some reshuffling of summation indices, we 
obtain 

Af _ f ; 

_q \ ^ jbc \ v ^aibd pa-at+d+c i » O o\ 

2 e-x iN a (.h-,£) 1^ J h+i ( A - 23 ) 

b,c,d£X i=l 

as the contribution to (|A.22[) which is proportional to a given er^ ^ . The multi- index 
a - a t + d + c= (ai, . . . , Oj_i,d, Oj+i, . . . ,o^,c) (A.24) 

in (j A.23|) is obtained from a by replacing the ith element by d and then appending c. 
This completes the proof. 
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